nb.MI <- function(d,outcome,add.pred.MI,add.pred,wald)	{

	wrap <- function(x)	{ paste0("$(", formatC(x,format = "f",3), ")$") }

	out.coef	<- vector(mode = "list", length = 10)
	out.var		<- vector(mode = "list", length = 10)
	out.sims	<- vector(mode = "list", length = 10)
	out.n		<- vector(mode = "list", length = 10)

	v <- 1
	for(v in 1:10)	{

		d.temp <- d

		##  expand predictors that are multiply imputed
		add.pred.MI.expanded <- ""
		m <- 1
		for(m in 1:length(add.pred.MI))	{
			add.pred.MI.expanded <- paste0(add.pred.MI.expanded, " + ", gsub("@",v,add.pred.MI[m]))
										}

		f <- as.formula	(
		paste0(	outcome, " ~ ", add.pred.MI.expanded,
				" + FirstSessionYear + Denunciation + Euthanasia + FinalPhase +",
				"ExterminationEinsatzgruppen + ExterminationCamps + ExterminationOther +",
				"Judicial + NSDetainment + NSOther + WarCrimes",
				add.pred)
						)


		##  fit model twice to find out whether any obs/covars
		##  need to be dropped because of missingness
		glm.out <- glm(f, data = d.temp, family = poisson(link = "log"))
		if(!is.null(na.action(glm.out)))	{ d.temp <- d.temp[-na.action(glm.out),] }

		hurdle.out 		<- hurdle(f, data = d.temp, dist = "negbin", link = "probit")
		out.coef[[v]]	<- coef(hurdle.out)[!is.na(coef(hurdle.out))]
		out.var[[v]]	<- vcovCL(hurdle.out, fix = TRUE, cluster = cbind(d.temp$Court,d.temp$FirstSessionYear))
		##  clustered at level of assignment

		##  consistency check
		stopifnot(length(out.coef[[v]]) == nrow(out.var[[v]]))

		##  record n
		out.n[[v]] <- nrow(d.temp)

					}

	##  consistency checks
	stopifnot(sum(is.na(out.coef)) == 0)
	stopifnot(sum(is.na(out.var)) == 0)
	stopifnot(length(unique(sapply(out.coef,length))) == 1)


	##  compute correct MI coefs and standard errors
	temp <- data.frame(matrix(NA,length(out.coef[[1]]),2))
	rownames(temp) <- names(out.coef[[1]])
	colnames(temp) <- c("coef","se")

	m <- 1
	for(m in 1:nrow(temp))	{

		out <- 
		pool.scalar(Q = unlist(lapply(out.coef, function(x)x[m])),
					U = unlist(lapply(out.var, function(x)x[m,m])),
					n = min(unlist(out.n)))

		temp[m,1:2] <- c(out$qbar, sqrt(out$t))

							}


	##  format results
	temp$z <- abs(temp[,1] / temp[,2])

	asterisks <- rep("", nrow(temp))
	asterisks[temp$z >= qnorm(0.975)]  <- "^{*}"
	asterisks[temp$z >= qnorm(0.995)]  <- "^{**}"
	asterisks[temp$z >= qnorm(0.9995)] <- "^{***}"

	temp[,1] <- paste0("$", formatC(round(temp[,1],3), format = "f", 3), asterisks, "$")
	temp[,2] <- wrap(round(temp[,2],3))


	##  Wald test
	Wald.out <- rep(NA,length(wald))

	if(!is.null(wald))	{

		##  judges
		Q <- matrix(0,1,length(out.coef[[1]]))
		Q[1,wald[[1]][1]+1] <- 1
		Q[1,wald[[1]][2]+1] <- -1
		Wald.out[1] <- unlist(MIwaldtest(qhat = out.coef, u = out.var, Cdes = Q, rdes = 0)$stat[4])

		##  prosecutors
		Q <- matrix(0,1,length(out.coef[[1]]))
		Q[1,wald[[2]][1]+1] <- 1
		Q[1,wald[[2]][2]+1] <- -1
		Wald.out[2] <- unlist(MIwaldtest(qhat = out.coef, u = out.var, Cdes = Q, rdes = 0)$stat[4])

						}

	return(list(temp,formatC(round(Wald.out,3), format = "f", 3), unlist(out.n)))

}







#

